Influence of vibrational modes on the quantum transport through a nano-device 
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We use the recently proposed scattering states numerical renormalization group (SNRG) approach 
to calculate I(V) and the differential conductance through a single molecular level coupled to a 
local molecular phonon. We also discuss the equilibrium physics of the model and demonstrate 
that the low-energy Hamiltonian is given by an effective interacting resonant level model. From 
the NRG level flow, we directly extract the effective charge transfer scale F ff and the dynamically 
induced capacitive coupling U c s between the molecular level and the lead electrons which turns 
out to be proportional to the polaronic energy shift E v for the regimes investigated here. The 
equilibrium spectral functions for the different parameter regimes are discussed. The additional 
phonon peaks at multiples of the phonon frequency uq correspond to additional maxima in the 
differential conductance. Non-equilibrium effects, however, lead to significant deviations between 
a symmetric junction and a junction in the tunnel regime. The suppression of the current for 
asymmetric junctions with increasing electron-phonon coupling, the hallmark of the Franck-Condon 
blockade, is discussed with a simple framework of a combination of (i) polaronic level shifts and (ii) 
the effective charge transfer scale r e ff. 

PACS numbers: 03.65.Yz, 73.21. La, 73.63.Kv 



I. INTRODUCTION 

In the quest for size-reduced and possible low-power 
consuming electronic devices the proposal 1 of using 
molecular junctions for electronics has sparked a large in- 
terest in understanding the influence of molecular vibra- 
tional modes onto the electron charge transfer through 
a molecule. In the simplest building block of molecular 
electronics, a molecule is connected to two leads. The 
non-linear current through such a device can be con- 
trolled by an external gate tuning the molecular levels^ 
In some cases a sudden drop of the current has been ob- 
served with increasing bias voltage- which translates into 
a negative differential conductance. Also, hysteretic be- 
havior of the I(V) curved has been reported when sweep- 
ing the voltage with a very small but finite rate. It 
has been suggested that such a reduction of conduc- 
tance and the hysteretic behaviour might originate in 
conformational changes in these complex molecules £ Vi- 
brational coupling has also been found of importance in 
break junctions 5 and suspended carbon nanotube quan- 
tum dots£~— An excellent reviewiS by Galperin et al. 
summarizes comprehensively the different theoretical ap- 
proaches and experimental findings. 

The theoretical description of such molecular junc- 
tions focuses only on those molecular levels and vibra- 
tionals modes which are relevant for the transport. In 
its simplest version 1 ^— a single level coupled to a lo- 
cal Holstein phonon has been considered. Typically rate 
equations 1 ^ or Keldysh-Green function approache a 10 : 14 
have been applied to this problem. One either expand 
the self-energy in powers of the electron-phonon coupling 
in the weak-electron phonon coupling regime 1 ^ or start 
from the exact solution of the local problem using the 
Lang-Firsov transformation 1 ^ and expand in powers of 
the tunneling matrix element^ The polaron formation 



on a molecular wire as a mechanism for negative differen- 
tial resistance 6 - - — which has been proposed uses a simple 
mean-field approximation^ Using the imaginary-time 
formalism ) 16 ' 17 however, the I-V curve of single orbital 
molecular junction does not show phononic site peaks j 1 ^ 
Whether this result prevails when the fit-function 1 ^ for 
the electronic self-energy is replaced by the exact solution 
remains an open question. Recently, the iterative path- 
integral approach has also been successfully applied^ to 
calculate quantum transport for moderate and high tem- 
peratures compared to the charge-transfer rate Tq. 

In this paper, we will briefly review the known physics 
of such a polaronic model from a renormalization group 
perspective. The low-energy Hamiltonian of the min- 
imal model for molecular devices is given by a effec- 
tive interacting-resonant level modeli^Sr— a considerably 
large Coulomb repulsion U e g is dynamically generated 
which governs the zero-bias transport as function of the 
gate voltage 2 ^ as well as the shape of the spectral func- 
tions, as we will demonstrate in our paper. We will 
present a detailed scaling analysis of the renormalized 
charge-transfer rate r e ff and U e g in the weak to inter- 
mediate coupling regime which covers a complementary 
regime of the once studied recently by Eidelstein et al^ 
We extend the discussion of the equilibrium properties to 
the single-particle spectral functions which governs the 
transport in the tunneling regime. Using the scattering 
states numerical renormalization group (SNRG)22r— we 
present the non-perturbative results for I-V characteris- 
tics of the model far from equilibrium at low temperature 
augmenting recent studies of other non-perturbative ap- 
proaches to larger temperatures^ 

Focusing on a single vibrational mode ) 10 ' 13 ' 19 ' 22 : 26 the 
equilibrium physics of two extreme limits have been well 
understood. 

In the adiabatic limit, where the phonon frequency is 
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the smallest energy scale of the problem a small electron- 
phonon coupling yields a reduction of the phonon fre- 
quency by particle-hole excitations. In leading order, 
the correction to the electronic self-energy is quadratic 
in the coupling constant. This limit has been pioneered 
by Caroli et al.— in the context of tunnel junctions and 
applied to molecular junctions using the self-consistent 
Born approximation^ 

In the opposite limit for very small tunneling rates t a , 
one starts from the exact solution of the local problem, 
t a = 0, by applying a Lang-Firsov transformation ! 1 ^ 28 
A displaced phonon with an unrenormalized phonon fre- 
quency cl> and a polaron with a shifted single-particle en- 
ergy is formed locally. When tunneling, the electron has 
to be extracted from the polaronic quasi-particle which 
can be done at many difference excitation energies differ- 
ing by multiples of the phonon frequency. Each pole of 
the single-particle Green function contributes only with 
a fractional weight to the spectral If the phonon en- 
ergy ujq is large compared to the electron energy scales, 
the bosonic mode can be considered in its ground state, 
which leads to an exponential renormalization of the 
tunneling matrix element t a — > t a exp(— g 2 /2), where 
9 = A p ft,/wo and \ p h denotes the electron-phonon inter- 
action strength. In this anti-adiabatic limit, the strong 
electron-phonon coupling yields a polaronic shift of the 
single particle level and, depending on the bare parame- 
ters, a reduction of the tunneling rate. This leads to the 
Franck-Condon blockade in quantum transport ^ 10 ' 13 ! 19 

In the limit of vanishing charge transfer rates, i. e. 
(r /uj ), (T /X p h) -C 1, the I-V characteristic exhibits 
rather sharp steplike features which are a reminiscence 
of the exact solution of the local spectral functio n 15 i 28 
as recently obtained using kinetioSS, equation, equation 
of motion decoupling schemes 3 ^ or master equations^ 
or a Keldysh Green function 3 ^ approach. Within such 
an approach negative differential resistance is found in 
regimes which are complementary to our scattering states 
NRG approach used in this paper. Such steplike features 
with negative differential resistance have been observed 
in suspended carbon nanotube experiments 

The problem becomes more difficult in the crossover 
regime between the adiabatic and the non-adiabatic 
regime. Often simple rate equations or self-consistent 
Born approximation^ fail to describe the proper renor- 
malization of the parameters. Recently, Wilson's nu- 
merical renormalization group (NRG) approac h 32 i 33 has 
been adapted to this problem^ A thorough study22 
has demonstrated the power of this non-perturbative ap- 
proach to reveal the interplay between the different en- 
ergy scales of the problem in the crossover regime. An 
extended anti-adiabatic regime has been identified where 
the bare charge transfer rate Tq exceeds the phonon fre- 
quency but remains below the polaron shift E p = g 2 ujQ. 
The scaling behaviour of renormalized charge transfer 
rate r e g- was obtained as function of the g and the po- 
laron shift E p . 

Galperin et alJ^ discuss the possibility of polaron for- 



mation on a molecular wire as a mechanism for negative 
differential resistance (NDR) uses a simple mean-field ap- 
proximation. Such hysteresis arise for branch cuts of the 
the non-linear equations at sufficiently larger coupling. 
It remains unclear whether this effect survives in the ex- 
act solution of the problem. A prominent example of a 
false hysteresis in the I-V characteristics is the conserv- 
ing GW-approximation to Anderson model out of equilib- 
rium as shown in a paper by Spataru et all 3 - 5 - In a recent 
paper it has been proposed 3 ^ that bistability signatures 
in non-equilibrium charge transport through molecular 
quantum dots might be linked to subtle differences in 
the initial conditions. 



A. Plan of the Paper 

The paper is organized as follows. We start with a 
definition of the model considered in this paper and re- 
late our form of the electron-phonon coupling to choices 
in the literature in Sec. Ill Al Decoupling the local level 
from the two leads, the local dynamics can be solved 
exactly by the well-known Lang-Firsov transformation 15 
summarized in Sec. Ill Bl We briefly review the numeri- 
cal renormalization group (NRG) approach to quantum 
impurities and the scattering states NRG used for the 
quantum-transport calculations in the Sees. Ill Cl and lll Dl 

The Sec.|nT]is devoted to present the equilibrium NRG 
results. We extract the parameters of the effective low- 
energy Hamiltonian in the regime ujq > T c g which is 
given by an interacting resonant level model (IRLM) . Af- 
ter summarizing the technical details how to extract the 
effective parameter directly from the NRG level flow, we 
discuss the scaling properties of r c ff and U c g as function 
of the bar model parameters in Sec. IHI CI 

The equilibrium spectral function already reveals im- 
portant information for the quantum-transport since it 
is proportional to the transfer-matrix in the tunneling 
regime. Therefore, we present equilibrium spectral func- 
tions for symmetric and asymmetric junctions as well as 
benchmarks for our non-equilibrium Green function al- 
gorithm in Sec. IIVI 

Scc.[V]is devoted to the results for the non-equilibrium 
quantum transport. We present our data for symmetric 
and asymmetric junctions for different phonon frequen- 
cies and junction asymmetries. We comment on the ob- 
servation of the Franck-Condon blockade physics and the 
recovery of the tunneling limit. We conclude with a short 
discussion and an outlook in Sec. I VII 

II. THE ANDERSON-HOLSTEIN MODEL 

A. Definition of the model 

The minimal, non-trivial model for molecular- 
electronic o 10 i 19 comprises a single spinless level whose 
charge is coupled locally to a single Holstein phonon 
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FIG. 1. (Color online.) Minimal model of a single-molecular 
level with energy Ed coupled to two leads with tunneling 
matrix elements th and tn. Depending on the local charge 
configuration jO) or jl), the Holstein phonon ground state is 
shifted. The relative displacement between the two configu- 
ration is given by the dimensionless electron phonon coupling 
g = Apfc/wo- The phononic excitations for a fixed charge are 
multiples of the oscillator energy luq. 



stemming from the dominating vibrational mode of the 
molecule. For the current transport this level is coupled 
to the two leads which are often considered as featureless 
bands for simplicity. In real materials, band features are 
important but only influence the single-particle proper- 
ties which can be accounted for in a frequency dependent 
charge transfer rate T(lj). 

This spinless Anderson-Holstein model is depicted 
schematically in Fig. [TJ Depending on the local charge 
configuration, the local harmonic oscillator is displaced 
and the distance between the two ground states is given 
by g (see below.) At the particle-hole symmetric point 
and in absence of a coupling to the leads, the displaced 
oscillator ground states are given by the two coherent 
states | ± g/2), and we immediately can understand the 
underlying Franck-Condon physics by the suppression of 
the overlap (—g/2\g/2) = cxp(— g 2 /2) of the two ground 
states with increasing electron-phonon coupling g. 

The spin-less two lead resonant level model (RLM) de- 
fined by the Hamiltonian 

Hq = H imp + Ht + Hi ea d s (la) 
H? mp = E d dU (lb) 

r t= £ -t=£(^+ c L<*) (ic) 

a—L,R v k 

Hi e ads = £ ^2 £kaC ka C ka (Id) 
a—L,R k 

is often use d 10 i 19 as the simplest model to describe quan- 
tum transport through a molecule or a suspended nano- 
bridge-* 3 -^ d(d') annihilates(creates) an electron on the 
device with energy Ed, and c ka creates an electron in the 



lead a with energy Ska- The local charge-transfer rate 
to each lead a is given by r a = nv^p a (0), where p a (u>) 
is the density of states of lead a. Throughout the paper 
we will use a lead independent constant density of states 
p a (uj) — (1/2D)Q(D — for simplicity; D denotes the 
band width of both leads. 

We only account for a single local vibrational mode 
with energy uo created by whose dimensionless dis- 
placement operator x = (tf + b) is coupled to the density 
fid = <$d of the local level via 

H ph = u tfb + \ ph (tf + b) (n d - M . (2) 

The form of the interaction in @ ensures particle-hole 
(PH) symmetry for E d — and particle-hole symmetric 
leads. 

Focusing only on the local dynamics defined by Hi mp — 
Himp + Hph , the harmonic oscillator will be shifted upon 
changing of the local occupancy away from half filling. 
This can be made explicit by introducing the arbitrary 
constant no and making the trivial substitution 

S d - 1 = (dU-n ) + (n Q - I) . (3) 

Then Hi mp takes the exact form 

H lmp = uo&b + X p h(b^ + b) (d)d - n ) 

+E d dU + 2E p (n ~^) 2 (4) 

where the displaced phonon is created by defined as 
&t = &t +5 ( no _I), ( 5) 

and g = \ p h/oJo- The new single particle energy Ed 

E d = E d + 2^(1 -n ) (6) 
luo 2 

A 2 

contains the polaronic energy shift E p = = u>og 2 - 
This polaron shift E p plays an important role of defining 
the different regimes^ of the model in addition to the 
dimensionless coupling constant g and the charge transfer 
rates r a . 

In a mean-field decoupling of the electron-phonon in- 
teraction, no would be replaced by self-consistent local 
charge expectation rid = (d^d) in Eqs. (j4][6]). A positive 
(negative) Ed leads to a depletion (filling up) of the lo- 
cal level where the effective level Ed is further shifted to 
higher (lower) energies by a term proportional to E p . 

This has a profound impact on the zero-bias conduc- 
tance as function of the detuning of the level Ed by an 
external gate voltage. The local average occupation nd 
and therefore, the total displace charge A./V due the pres- 
ence of the impurit y 38 ' 39 is initially controlled by the ratio 
Ed/T e ff, where T g is the low temperature charge fluctu- 
ation scale in the Fermi- liquid fixed point. In the strong 
coupling regime, g ^ 1, the effective level Ed becomes 
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almost independent of the gate voltage controlling the 
barelevel Ed- Hence, the occupation will be nearly inde- 
pendent of Ed and, therefore, the conductance will have 
a plateau unless \Ed\ exceeds the polaron shift E p . The 
zero-bias conductance depicted in Fig. 13 of Ref. [22| is a 
direct consequence of Eq. (|6|) . 

We also can use the parametrization of Hi mp with tiq as 
defined in Eq. §4§ to connect Hi mp with the local molec- 
ular Hamiltonian H m 

H m = uotfb + \ ph (tf + b)d ] d + edU (7) 

often used in the literature! 10 ' 11 ' 19 Setting n = 0, e = 
Ed + E p and neglecting the ground state energy shift, 
Himp becomes identical to H m . 

Keeping e fixed in H m and increasing the electron- 
phonon coupling X p h as in Ref. [l9| translates into a renor- 
malization of the single-particle energy by a polaron shift 
Ed = s — E p in Hi mp . Starting from e = 0, this implies 
a population increase of the fermion level upon increase 
of Xph . In this parametrization it becomes apparent that 
Xph induces a detuning away from particle-hole symme- 
try. 

Combining the RLM as given in Eq. (TTJ with a coupling 
to the local Holstein phonon to H — Hq + H p h defines 
the spinless Anderson-Holstein model. 



given by 

H' = H' lmp + J2 £kaC L Ck <* 

a=L,R k 

+E^E( e9(b " b) ^+4>- s(bt - b) ) 

a V iV k 

Hl mp = E d tfd + u> b% (10) 

after the unitary transformation H' = U^ f HUlf. We 
have dropped the bar from the transformed phonon b = 
U^pbULF in Eq. (flUl) to keep the notation simple. While 
the local Hamiltonian H' imp became simple and diagonal, 
we acquired a complicated tunneling term with an expo- 
nential electron-phonon coupling e ±9 ^ b ~ b K 

Although we have always used the original Hamilto- 
nian in our NRG calculations, the transformed Hamilto- 
nian of Eq. (|10|) is very convenient to gain some deeper 
insight into the effective low-energy Hamiltonian gener- 
ated by the renormalization group transformations. This 
effective model is discussed in detail in Sec. IIIIAl below. 



B. Lang-Firsov transformation 

This Anderson-Holstein model is well studie d 10 ' 19 ' 22 ' 34 
and a text-book example for an exactly solvable model 28 
in the limit t a = 0. Already in the 1960s, it was shown 
that the local Hamiltonian H imp — H® mp +H p h can be ex- 
actly diagonalized using a Lang-Firsov transformation, 15 

U Lp = e -9V-b)tfd-i) (g) 

The new elementary excitations are polarons annihilated 
by the operator d 

d=Ul F dU LF = e- g ^- b] d (9) 

and a free phonon with unrenormalized energy ujq. 
The polaron shift E p enters the ground state energy. 
Acting on the bosonic vacuum, the operator Ulf — 
e -s0 1 -&)(<i 1 <2-±) g enera tes a coherent state |±g/2) which 
describes the ground state of a displaced harmonic os- 
cillator by the dimensionless displacement Aa; = ig/2 
depending on whether a local charge was present or ab- 
sent (see also the schematic figure [TJ) This reflects the 
underlying Franck-Condon physics extensively discussed 
in the literature) 9 - ' 10 ' 13 

Neglecting the ground state energy shift, the total 
Hamiltonian H' of the impurity coupled to the leads is 



C. Numerical renormalization group and the 
Anderson Holstein model 



The properties of quantum impurity systems such as 
defined by H = Hi„ lp + Hp + Hi ea d s can be very ac- 
curately calculated using the numerical renormalization 
group. At the heart of this approach is a logarithmic dis- 
cretization of the continuous bath, controlled by the dis- 
cretization parameter A > L 32 i 33 The continuum limit is 
recovered for A — >■ 1. Using an appropriate unitary trans- 
formation^ the Hamiltonian is mapped onto a semi- 
infinite chain, with the impurity coupled to the first chain 
site. The Ath link along the chain represents an expo- 
nentially decreasing energy scale: ~ K~ N / 2 . Using 
this hierarchy of scales, the sequence of finite-size Hamil- 
tonians V. n for the A-site chain is solved iteratively, dis- 
carding the high-energy states at the conclusion of each 
step to maintain a manageable number of states. The 
reduced basis set of T-Ln so obtained is expected to faith- 
fully describe the spectrum of the full Hamiltonian on a 
scale of Dn, corresponding to the temperature Xjv ~ fljv- 
Details can be found on the review^ 3 - by Bulla et al.. 

Hewson and Meyer— pioneered the application of the 
NRG on the single-lead version of Hamiltonian H = 
H + Hph- The standard NRG discretisation^ 3 , of 
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Hieads maps the model onto a chain Hamiltonian 
H N A- {N - 1)/2 = E d Sd + u) tfb + Aph(6 f + b) (did - 



N 



enters as the expansion coefficient where 5 > is an 
infinitesimally small energy scale required to select the 
correct boundary conditions. 

Defining t = y/t 2 L + t 2 R , the ratio r R[L) = t R[L) /i 
denotes the relative tunneling strength of each lead a 
to the impurity d- level. The resonant level self energy 
A a (uj + iS) in Gq(cj) is given by 



WA-" /2 e f f f 

/ j / j 11 z-naJnaJn 



(11) 



a n— 
N-l 



a n— 

using a NRG parameter A > 1. t n ,e n = 0(1) and e„ = 
for PH-symmetric leads. The dimensionless parameters 
E d ,uj , X p h are related to original parameter a 32 ' 33 via the 
scaling factor s = D(l + l/A)/2: E d = E d /s,u} — 
loq/s, Xph — Xph I ' s and the renormalized tunneling V a — 
t a /s. Using a suitable number of bosonic excitions Nb, 
Hn is then iteratively diagonalized. A much more de- 
tailed introduction to the NRG can be found in the Re- 
view by Bulla et 



D. Scattering-states numerical renormalization 
group approach to quantum transport 

For the calculating of the current /(V) through 
the molecular level at finite bias V across the two 
leads, we have used the recently proposed scattering- 
states numerical renormalization group 2 ^— (SNRG) ap- 
proach to quantum transport. The SNRG is based 
on an extension of Wilson's NRG to non-equilibrium 
dynamics.— ~— Using the time-evolved density operator, 
the non-equilibrium steady state retarded Green func- 
tion can be calculated^ 3 - Below, we only give a rather 
brief summary of the approach. A more detailed deriva- 
tion and discussion of the method can be found in Refs. 

MM. 



1. Definition of the scattering states 

In the absence of the electron-phonon interaction the 
RLM defined in Jl} can be diagonalized exactly in the 
continuum limi t 16 ' 44 " — by the following scattering-states 
creation operator o 23 ' 24 ' 47 ' 48 



7 t 

lea. 



+ t a y/p a (e)Go(e + iS) 

jtg' y/Pa 1 (e') f 



(12) 



dt + £ / d£ ' 



e + i5 — e' 



a = L(R) labels left (right) moving scattering states cre- 
ated by l\,y L ^ ■ In this equation, the local retarded 
resonant-level Green function of the d-level 



G5H 



+ iS - E d - A(lo + iS)}" 1 



(13) 



A(w + iS) 



de- 



iS 



A ^ tn(fn a fn+la + /n+lce/na) 



»e[A(w)] - ir(w) 



(14) 



and its imaginary part T(u) denotes the energy depen- 
dent charge-fluctuation scale which becomes a constant 
in the wide-band limes of a featureless band. 

In the limit of infinitely large leads the single-particle 
spectrum remains unaltered, and these scattering states 
diagonalize the Hamiltonian (flj 



Ho = H(Xp h = 0) 



E 

a=L,R 



(15) 



up to a neglected ground state energy shifts 

The creation operator jj a is a solution of the opera- 
tor Lippmann-Schwinger equation 4 ^ and, therefore, the 
corresponding state break time-reversal symmetry. The 
necessary boundary condition for describing a current- 
carrying open quantum system is encoded in the small 
imaginary part +iS entering Eq. (|12M14p required for con- 
vergence when performing the continuum limit Vol. — > 
oo. Since left and right movers at the same energy are 
time-reversal pairs, time-reversal symmetry is restored at 
zero bias yielding a exactly vanishing net current in that 
limit. 

To avoid possible bound states, we will only consider 
the wide-band limit: D ^> max{|e<j|, r Q , |V|}. Further- 
more, r = r(0) = Tl + Tr ; is used as energy unit in this 
paper. We measure the coupling asymmetry by the ratio 
R = Tt,/Tr. A perfect unitary limit of e 2 /h conductance 
quantum can only be reached for R — 1 with R — > oo and 
R — > correspond to the tunneling regime. 

Hershfield has shown that the steady-state density op- 
erator for a current-carrying non-interacting quantum 
system retains its Boltzmannian form 4 ^ 



-P(Ho-Y ) 



Po 



Tr e -/3(«o-Vb) 



Yo = ^2po 

OL<J 



"fso-aYscra 



(16) 



at finite bias. The Y$ operator accounts for the different 
occupation of the left-moving and right-moving scatter- 
ing states. p a denote the different chemical potentials of 
the leads. Since H is bilinear, the transport is perfectly 
ballistic and the total current is given by the difference 
between the current to the left and the current to the 
right. 
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All steady-state expectation values of any operator are 
calculated using p for the non-interacting problem which 
includes the finite bias. In the absence of the electron- 
phonon interaction this is a trivial and well-understood 
problem. It was shown 4 ^ that the current obtained with 
this density-operator po is identical to the current cal- 
culated using a generalized Landauer formula based on 
Keldysh Green functions <£2r— 

This form of po stated in Eq. (| 16[) remains valid even for 
the fully interacting system^ 5 - when replacing Hq — s- H, 
and replacing Yq Y. Since the Y must be constructed 
from the many-body scattering states, its explicit ana- 
lytical expression is unknown for a general Hamiltonian 
H of an interacting systems. 

We, therefore, proceed in two steps using the argu- 
ments outlined in the literature At first, we add to 
H a fictitious electron-phonon interaction term H' e _ ph 
which commutes with Yq, i.e. [H' e _ ph ,Y ] = 0. Hersh- 
field's argument also yields a steady-state density opera- 
tor of the form p' = exp[-^(K + H' e - P h ~ %)\l z - For 
the model under consideration we have chosen 

H' e _ ph = w tfb + \ ph (tf + b)l J2 r l d l d « - \ (17) 

\a=L,R J 

where the d a are the operators obtained from inverting 
Eq. CE} 



d a = t de^p(e)[Gl{e)y le 



(18) 



which fulfill the anti-commutation relation {d a ,d^} = 
8 a p. The annihilation operator d of a local electron 
on the device is reconstructed by the linear combina- 
tion of its left-mover and right-mover contributions d = 
TRdn -Vrj^dj,. We note that H' e _ ph approaches H e ^ p h in 
the extreme tunneling limit of R — > oo (or R — > 0.) 

In a second step, we perform the time evolution of 
p' with respect to the fully interacting Hamiltonian to 
infinitely long time: the density operator p(t) progresses 
from its initial value p' at t = as 



(19) 



where we set h—1. 



2. Scattering- states numerical renormalization group 

The basic idea of the scattering-states numerical renor- 
malization group (SNRG) approach is summarized as fol- 
lows. 

(I) Knowing the analytical form of the non-equilibrium 
density operator p' , we can discretize scattering states 
on a logarithmic energy mesh identically as in the stan- 
dard NR G 23 i 33 and perform a standard NRG using 
Kq = V.(X p h = 0) + H' e _ ph — Y a . The density operator 
p' (V) contains all information about the current carrying 



steady-state for the Hamiltonian Kq + F 4^ 

(II) Starting at time t = 0, we let the system evolve 
with respect to the full Hamiltonian Hf = H{\ p h > 0) 
Then, the density operator p{t) progresses from its initial 
value /5q at t = according to Eq. (TTO)) . Since we quench 
the system only locally, it is a fair assumption that p(t) 
reaches a steady-state at t — > oo independent of initial 
condition for an infinitely large system, since all bath 
correlation functions must decay for infinitely long times. 

The finite size oscillations always present in the NRG 
calculation— are projected out by defining the time- 
averaged density operator 



lim — 

r— >oo T 



dtp(t) . 



(20) 



As a consequence, only the matrix elements diagonal in 
energy contribute for t — ¥ oo in accordance with the 
steady-state condition 



[n f ,p 00 ] = o 



(21) 



Even though poo remains unknown analytically, we 
can explicitly construct it numerically using the time- 
dependent NRGia-i2 (TD-NRG.) 

(Ill) The steady-state retarded Green function is de- 
fined as 



G r AiB (t) = -m Poo [A(t),B} s e(t), 



(22) 



where A(t) = e m f t Ae~ iH f t 1 [A(t),B] s denotes the com- 
mutator (s = — 1) for bosonic, and the anti-commutator 
(s = 1) for fermionic correlation functions. This Green 
function can be calculated using the time-dependent 
NRG 4 -4i and extending ideas developed for equilibrium 
Green functions^ The details of the algorithm is derived 
in Ref . IS 



3. Current as function of the bias voltage 

The current I a is defined as a charge current^ from 
the lead a to the local d-level. It has been show n 45 i 47 i 51 
that for the model investigated here, the symmetrized 
current 



I = rllL ~ r 2 L I R 



(23) 



is related to the steady-state spectral function of the local 
level, pd(ui,V) = Sm[G£(o; — i0 + , V)]/w by a generalized 
Landauer formula 

I(V) = -± / du \f R (u) - f L {u)] T 7rp d (u, V)(24) 
where f a (oj) = f(uj — p a ). The prefactor 

e 2 4r L r fl . 



Gn = 



h Tl 



(25) 
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measures the asymmetry of the junction. Go reaches 
the universal conductance quantum e 2 /h for a symmet- 
ric point-contact junction, i. e. Tl = Tr, and is strongly 
suppressed in the tunneling regime T a <C r_ Q . The 
two chemical potentials [x a are set to /il = —r R V and 
)ir = r\ V as function of the external source-drain volt- 
age V consistent with a serial-resistor model which is 
required by the current conservation Ir = —II- 

While the conserving Keldysh Green function ap- 
proach treats the problem of initially decoupled leads 
and propagate the system to a current-carrying steady 
state, the scattering states NRG starts from an initial 
current carrying steady state of the non-interacting prob- 
lem. Employing the TD-NRG 40 ^ we let the system 
evolve by calculating the full density operator p(t) nu- 
merically in the limit of t — > oo using the Hamiltonian of 
the interacting system. 

In a continuum limit, any correlation function will 
have a finite correlation time above which the memory 
of this initial state will be lost.— ^ As a consequence, 
a unique steady state is approached which is compati- 
ble with the imposed boundary condition as pointed out 
by Hersfieldi 45 ' 53 Our starting point is therefore identi- 
cal to any approach using a propagation along a Keldysh 
contour^ The differences arise from the inevitable ap- 
proximations made in Keldysh perturbation theory by 
selecting a subclass of diagrams for a practical calcu- 
lation of the non-equilibrium self-energies. In contrast, 
our approach does not make any of such approximations 
and, therefore, includes the full interaction to infinitely 
high order. It only comprise a systematic but well con- 
trolled error stemming from discretization of the bath 
continuum^ii 2 .^ inherent any NRG (or DMRG^) ap- 
proach. 

A hysteretic behavior 4 is observed in some molecular 
junctions when sweeping the voltage a very slow but fi- 
nite rate. Using a simple mean-field approximation, the 
possibility of polaron formation on a molecular wire has 
been proposed as a mechanism for the observe NDR. 11 
Such hysteretic behavior in self-consistent equation can 
also occur in the employed approximations when increas- 
ing the coupling constant beyond the validity range of 
such approximations. Such many-valued solution might 
not survive in an exact solution of the problem and their 
robustness need to be checked against the inclusion of 
higher order contributions. A prominent example of such 
a known false hysteresis in the I-V characteristics is the 
GW-approach to Anderson model out of equilibrium^ 

The theoretical description of such hysteretic behavior 
requires tracing the same experimental conditions in the 
simulations. Starting from the initial condition of a fully 
interacting current-carrying steady state for a given bias 
voltage the time dependent current I(V(t),t) must be 
calculated for a given rate of voltage change. 

In a recent paper by Muehlbacher et al^ it has been 
proposed that signatures of bistability in non-equilibrium 
charge transport through molecular quantum dots might 
be linked to subtle differences in the initial conditions. 



This, however, is clearly beyond our scattering-states ap- 
proach which targets exclusively the stationary steady- 
state limit. 



III. EQUILIBRIUM RENORMALIZATION 
GROUP APPROACH 

Before we present the results for the quantum trans- 
port in Sec. |Vl we briefly review the equilibrium param- 
eter flow of the model. This provides the necessary un- 
derstanding of the energy scales and defines the different 
regime which become relevant for the quantum transport. 



A. Low energy Hamiltonian 

While Hewson and Meyer investigated the equilibrium 
dynamics of the more complex spin-full model^ 4 - using 
the NRG, we focus on the simpler spin- less model in this 
work. Its low energy fixed point is also given by an ef- 
fective resonant level model. In the limit loq > D,To, 
its effective tunneling matrix elements i° ff are estimated 
by using H' defined in Eqn. (flU)) : the local phonons are 
approximately in their ground state, and 

t a -> tf « t a (0\ exp(-«?(6 t - 6))|0) 

= ta eM~) • (26) 

In leading order, we obtain the renormalized charge fluc- 
tuation scale T c g roexp(— g 2 ) in this limit. 
By expanding the factor 

e -ff(6+-6) = i - g (bt - b) + O (g 2 ) (27) 

for small g <C 1 added by the Lang-Firsov transformation 
to the tunneling term 

H' T = t a (e°^-»tfc 0a + clJe-^-V) , (28) 

a 

it is apparent that an effective repulsive Coulomb inter- 
action term 

Hu=J2u e s(dU-^)(n 0a -^) (29) 

a 

is obtained in second-order perturbation theory in g. Of 
course, this type of effective interaction and all possible 
higher order contribution will be automatically generated 
by RG transformation in each NRG iteration step. Tak- 
ing into account the renormalization of the tunneling ma- 
trix element t a — > i° ff , we conjecture that U e g generated 
in the RG transformation is given by 




in leading order. Then, the dimensionless Coulomb re- 
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pulsion reads 



D 



Wo 



U a 



D 



-/(*) 



(31) 



where f(x) is an unknown scaling function of order 
O(l) which accounts for higher order corrections in the 
expansion (j2"7| . The analysis of the NRG data, pre- 
sented below suggests that the scaling variable is given 
by x — a(u>o/To)g 2 in the anti-adiabatic regime ujq 3> IV 
The function a(y) is a slow varying function of the order 
0(1) as shown in Sec HlTC] 

A perturbative treatment^ 2 - to second order in t a , pre- 
dicts U c s/D ps (4T()Ep)/(itujq) in the weak coupling 
limit, i. e. g -C 1, and U e g/D ps (4Tq/E p ) in the strong 
coupling limit, g ^> 1. Note that these results are only 
valid as long as Tq is the smallest energy scale and, there- 
fore, U e s/D remains a perturbative correction in the va- 
lidity range of the perturbation theory. This analytic re- 
sult agrees nicely with our conjecture of Eq. (f3~Tj) which 
will be backed up by our extensive NRG study below. 

Note, however, that the values for U e g/D extracted 
from the full NRG calculation can result in magnitudes 
of O(l) and, therefore, this regime would be beyond the 
reach of a second-order perturbation theory. Neverthe- 
less, the leading-order scaling with the model parameters 
remains well captured by Eq. (|3 1 1) as long as ujq/Tq > 1. 



B. Renormalized parameters 

In equilibrium, the anti-binding combination of lead 
operators, Co- = t^cqr — trCql decouples from the im- 
purity, and we are left with an effective one-band model 
since the local d-orbital is only connected to orthogonal 
combination cq + — trCor+tlCql- Therefore, we focus on 
an effective one-band model but with the charge-transfer 
rate T + T L + T R . 

Eidclstcin ct ali22 have already pointed out that the 
low-energy Hamiltonian of the spinless one lead Ander- 
son Holstein model can be mapped onto an effective in- 
teracting resonant level model (IRLM) 



^jIRLM _ jj 



Hi, 



imp 

A,- 

+*eff(d t C + c f d) 
H imp = ed<$d + U c g{d)d - ~)(cjc - i) 



(32) 



for effective band width D e g < luq. Its parameters Ed, t e s 
and U e g can be extracted in different ways from the NRG 
level spectrumj 22 i 32 ~ — 

In Ref. I22I . the effective parameter r o g oc i 2 ff is ex- 
tracted from the charge susceptibilty at the particle-hole 
symmetric point. We, however, use a different approach 
which allows for a definition of T e s for arbitrary values 
of Ed, away from the particle-hole symmetric point by 



applying the procedure proposed by Hewson et al.<££ 

In this method the renormalized parameters are ex- 
tracted directly from the NRG level flow. It is based 
on the analytically known low-energy stable Fermi-liquid 
fixed points and the exact solution of the RLM for a 
given NRG chain. The discretized effective RLM for the 
fixed-point dynamics, 



H RLM A -(N-l)/2 = gddU + + f t d) 

N 

+ £A-"/ 2 £ -„/t/„ (33) 



n=0 
N—l 



^A-«/ 2 f„(/t/„ +1 + /t +1 / n ) 



n=0 



is exactly diagonalize d 32 ' 33 to 



2 

h rlm = (eyelet + e^hjhi) (34) 



1=1 



for a NRG chain with odd number of lead sites - N even, 
et (hi) creates the Z-tli elementary particle (hole) excita- 
tion with the energy e Pj ; (sh,l-) Obviously, the impurity 
Green function (z) 



z-e d A( N -^/ 2 -[V eS } 2 A N ^g (z) 



(35) 



must have poles at these single-particle excitation ener- 
gies z = £ e ,i,£h,i- go(z) is the Green function at chain 
site n = in the absence of the impurity and can be 
calculated by a continuous fraction expansion^ 

The first single-particle excitation, Ei p , and the first 
single-hole excitation of the NRG level spectrum, E\h re- 
spectively, are very good approximation of s Pi i and en,i 
close to the Fermi-liquid fixed point. These two NRG 
energies are sufficient to determine the two unknown ef- 
fective parameters e d and V e g by solving the two coupled 
algebraic equations^ 



E^A-^- 1 ^ 2 e d (N) _ A^- 1 )/* 



2r eff (A) 2r eff (A) 7T 
-E^A-^-V/ 2 e d (N) A^- 1 )/ 2 



2r eff (A0 



2r eff (iV) 



g {E lp ) (36) 
go(~E lh ) 



for the poles of the Green function, where T c g = 
7rf 2 ff/ o(0). Note, that these are dimensionless parameters, 
all given in units of D(l + l/A)/2. 

In order to calculate U c s, we analyze the four eigen- 
states of iJQ RLM comprising only the local level and Jo- 
lt is easy to see that U e e is related to the energy differ- 
ence between the lowest particle-hole excitation and the 
sum of a single-particle and a single-hole excitation. The 
adaptation of equation (18) in Ref. [HI 

E lp + E lh - E lph = 2U cS (N)A^ N ~ 1 ^ 2 

x (hMl)| 2 Kfc(l)| 2 + l^,h(l)| 2 |^ c ,p(l)| 2 ) (37) 
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requires the knowledge of expansion coefficients ipd,pQ) 
of the annihilation operator d 



(N+2)/2 



d= ^ {^d,p{l)e p ,i + ^d,h{l)h\ ht 



1=1 



and ip c ,p(l) of the first lead chain site n = 

{N+2)/2 
1=1 



(38) 



(39) 



in terms of the elementary excitation of the low-energy 
fixed point Hamiltonian (|34l) . The coefficients \ipd. P (l)\ 2 
are determined by the weights of the poles of the finite 
size chain impurity Green functions 



(JV+2)/2 

E 

i=i 



|^, P (0I 2 , \i>dM 



z - e p j 



z - eh i 



(40) 



and coefficients \ipc,h{l)\ 2 by the corresponding Green 
functions of the first lead site n = 0, G^ (z), which is 
related to the local T- matrix V 2 G% (z) 



N-1/~iNi 



(41) 



= 9a(z) 



g«ff A (JV-l)/2 



Z - £ "f A(JV-l)/2 _ [V eS ]^A^g (z) 



Using the calculated parameters £d(N) and T (N), the 
spectral weights at the pole e p i are approximately given 
by 



iiMi)r = r 



1 



[V eS ] 2 A"-ig' (E lp ) 

E lp - EdA^-W 



^ )| a ^ ii#TO (42) 

and analog for Eh.i- We have replaced the true single- 
particle (single-hole) excitation energy e p .\ [sh,l) by the 
first NRG particle (hole) excitation energy E\ p (Eih.) 
Thereby, g' (z) denotes the derivative of go(z). 



C. NRG analysis of the renormalized parameter 



0.8 
0.6 
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— r e,f (N)/r 

-- U d '(N)/D 
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N 



FIG. 2. (Color online.) Example NRG flow for U(N)/D 
and T c g(N) vs NRG iteration N for the parameters u>o/To = 
Aph/Fo = 5. 
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FIG 3. (Color online.) (a) Anti-adiabatic and crossover 
regime: T e g as function of g 2 — (\ p h/uo) 2 for ujo/^o = 

2 2 

1,5,10, (b) rescaled Y e se 9 vs g ; (c) r e ff/ro vs E p in the 
adiabatic regime. The correction factor a(ujo = 0.2) = 
1.02, a(wo = 0.3) = 0.93. In the anti-adiabatic and the 
crossover regime, (a) and (b) N 3 = 1500, in the adiabatic 
regime (c) we kept N s = 2500 NRG states for accuracy. 



In this section, we present the results for the renor- 
malized parameter for the low-energy equilibrium Hamil- 
tonian extracted from the NRG level spectrum of the 
single-lead model as described in section IIII Bl We have 
used A = 1.5, kept N s = 1500 NRG states after each 
iteration, and included the lowest Nh = 400 local phonon 
states. The band-width of a lead with constant DOS 
p(u) = B(D - M)/2D was set to D = 100r , and T 
serves as energy scale throughout the paper. In this sec- 
tion, we only investigate the PH-symmetric limit by set- 
ting E d = 0. 

In Fig. H the flow of r cff (iV) and U cS (N) is depicted 
with respect to the NRG iteration N for one specific pa- 
rameter set, uj /T = Ap/j/To = 5. Fast convergence^ 
is achieved when approaching the low-energy fixed point 
for N — > oo. We have extracted the fixed point val- 
ues of £^ = limjv- 



ef(N), r cff = lim w ^ oo r cff (iV), 
U c s = liniAr^oo U e ff(N) for various values of the electron- 
phonon coupling \ p h and ujq. Since, e^ ff = for all PH- 
symmetric parameters, we do not show any results. In 
the ph-asymmetric regime, e^ ff ^ measures the effective 
PH-symmetry breaking external field. 
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The results for r e ff/r vs g 2 are shown in Fig. [3^a) 
for three different phonon frequencies uj /T = 1,5,10. 
For To <C cjq , we approach the anti-adiabatic limit 
in which the phonons can almost instantaneously re- 
act on the charge fluctuations. The data suggests that 
r c ff ~ To exp[— g 2 ] for loq — > oo and weak and moder- 
ate values of g. The closer ujq approaches the charge- 
fluctuation scale To, the stronger the deviations for this 
approximation. In the context of the spin-full model, 
Hewson et ali^i have already proposed a modification of 
the exponent exp(— g 2 ) — > cxp(— g h(\ p h, u) )) by a func- 
tion h which accounts for the additional correction. For 
large ujq, the leading contribution to h is a constant ap- 
proaching 1, however, we were not able to find a univer- 
sal scaling function h using our data which describes the 
crossover from weak to strong coupling, g < 1 and g > 1 
respectively . 

Restricting to the strong coupling limit, however, 
Eidelstein et ali^ 2 - were able to show that all graphs 
r e ff((?,u;o) collapse onto a single curve by introducing 
the universal function through the expression r c ff/r = 
exp[— g 2 F(R)]. The dimensionless parameter R = To/E p 
emphasizes the importance of the polaronic energy shift 
E p for the physics of the model. 

The validity range of such a universal function F(R) 
requires JJ « 1. For R — > 1, however, significant devi- 
ations from a universal scaling are already observed in 
Ref . U 

Since we focus on the regime R > 1, the study by Ei- 
delstein et al. 22 investigates the complementary strong- 
coupling regime to our weak and intermediate coupling 
regime. Although the equilibrium NRG can reach R <C 1 
as recently demonstrated, the TD-NRG and the SNRG 
are restricted to moderate values of g due to the increas- 
ing of discretization artefacts i 42 ' 56 

In the adiabatic regime, ujq <C To, field-theoretical 
arguments^- have been employed to show 

r -~ r °H§) < 43 » 

for max{E p /To, g, Ed/To} <C 1. Our numerical analy- 
sis for the adiabatic regime agrees with that perturba- 
tive result, but we find a pref actor of 1/2 up to very 
small corrections instead of 2/n stated in Eq. (|43j) . as 
shown in Fig. EJc). In the analytical calculation-^ the 
local Lorentzian spectral function was approximated by 
a simple constant p while in our numerics the full en- 
ergy dependency enters the calculations which probably 
accounts for the difference. Our data provides evidence 
for a crossover from a perturbative [1 — a rf ] correction 
factor to Tq with a ~ 1/2 in the adiabatic regime to an 
exponential suppression factor exp[— g 2 h(\ p h, ujq)] in the 
anti-adiabatic regime. 

Let us now discuss U e g dynamically generated by the 
electron-phonon interaction. The plot in Fig. [5] demon- 
strates clearly that such an effective Coulomb repulsion 
U(N) can be extracted directly from the NRG level flow 
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FIG. 4. (Color online.) (a) u>oU e s as function of g 2 for four 
different values of ojo/Tq = 1,5, 10,20. (b) f(x) — U c s/U an a 
as defined in Eq. (|31[) vs a(uio)g 2 , where we used the scal- 
ing factor a(w ) = 0.85,1,1.1 for uj /Y = 5,10,20. NRG 
parameters as in Fig. [3] 



via Eq. (|37[) and converges rather quickly to a low tem- 
perature fixed point value J7 e ff. 

The results for ujo U e g are depicted in Fig.QJa) for four 
different values of ujq- Qualitatively, the results can be 
understood by the conjectured form 

E4na _ 71" r e ff 2 _ n F e gE p 

D ~ 2~uj^ 9 ~ 2 w 2 

of Eq. ([3l"jl above where U c s cx J7 ana - u>oU e g increases 
linearly with x = g 2 and is exponentially suppressed for 
large x due to the renormalization of the bare hopping 
constant. Focusing only on the three graphs in the anti- 
adiabatic regime cu Q > r , we plot the ratio U e g/U ana as 
function of a(u>o)x. 

By an appropriate dimensionless scale a(u>o) which 
only depends on the phonon frequency, we are able to 
rescale x such that all ratios U e s/U an a collapse onto one 
universal scaling curve of the order O(l) in this regime 
as shown in Fig. |U[b). 

The simple arguments leading to the estimated analyt- 
ical value Uana do not hold in the extended anti-adiabatic 
and in the crossover regime, and this mapping onto the 
universal scaling curve fails for ojq/Tq = 1. As indicated 
already in Fig.[4ja), U e g is much larger and the additional 
renormalization for r c g- by the finite U c g in the effective 
IRLM might have to be taken into account when crossing 
over to the adiabatic regime wq < Tq. 



11 



IV. EQUILIBRIUM SPECTRAL FUNCTIONS 

A. Technical details 

Since the equilibrium spectral functions are used to 
calculate the zero-bias conductance within the linear re- 
sponse theory, we present some local spectra to set the 
stage for the non-equilibrium transport. Furthermore, 
we also can use the direct calculation of the equilibrium 
spectral function as a benchmark for testing the qual- 
ity of the non-equilibrium spectra obtained by the time- 
dependent NRG approac h 40 i 41 to non-equilibrium Green 
functions. 43 

For two trivial limits, exact solutions are analytically 
known, (i) In the absence of the coupling to the leads, 
t a = 0, the exact solution of the local spectral function 28 
comprises of a set of equidistant <5-peaks separated by the 
phonon- frequency ujq. Their spectral weights are given by 
modified Bessel functions whose arguments are tempera- 
ture dependent, 28 (ii) In the opposite limit, X p h = 0, the 
spectral function of is simply given by the Lorentzian of 
the RLM which will be weakly modified for < g <C 1. 

We have used the complete Fock space 
algorith m 43 i 54 ' 61 to calculate the NRG spectral function. 
The (5-functions of the spectral function have been 
replaced by the standard NRG broadening function 

%±K|)4^ffle-^*/-lf (44) 
b\aj n \y/iT 

characterized by the broadening parameter b 
and additionally averaged over N z different band 
discretizations! 40 ' 41 ' 62 The equation of motion^ 3 - exactly 
relates the self-energy 

EP " W " X »'M (45) 

to ratio of the correlation function 

F(z)=<(b + tf)d\eft2>(z) (46) 

and the local Green function Gd(z) =<^ji\(V^> (z). 

For Xph — 0, the correlation function F(z) vanishes: 
it is generated in first order by the electron phonon cou- 
pling. Therefore, F(z) cx X p h in leading order which im- 
mediately reproduces the leading order of magnitude of 
S ph (z) oc Xp h for the weak coupling regime X p h/To <C 1 
obtained from perturbation theory 2 ^ where the leading 
order Feynman diagram is depicted in Fig. [SJ 

The spectral functions discussed below are all obtained 
from G d {z) = [z - E d - A(z) - E ph (z)]- X where the 
total self-energy is given by the sum of S ph (z) and A(z) 
defined in Eq. (|14[). In our numerics we evaluate the 
Green function very close to the real axis and Pd{^>) — 
^mGdi^J ~ iS)/ir using a very small value for S/Tq — 

io- 10 -io- 7 . 




FIG. 5. (Color online.) Equilibrium spectral functions 
for Ed = 0, loo /To = 1 and various values of Aph/To = 
0.5, 1.5,2,3. NRG parameter: iV s = 1200, A = 1.5, N b = 80, 
N z = 512, 6 = 0.03. 



B. Particle- hole symmetry 

1. Crossover regime 

In Fig. [5j the positive part of the symmetric lo- 
cal spectral function pd(w) — 3mG,j(u) — i5)/t: is de- 
picted for a series of electron-phonon couplings X p h /Tq = 
0.5, 1.5, 2.0, 3.0 and a fixed phonon frequency loq/Tq = 1. 
For Xph/To = 0.5, the RLM Lorentzian is broadened 
and a kink at to w loq is related to a sharp rise of the 
self-energy at to. With increasing X p h, more and more 
phonon replica of the resonance peak to = Ed = are 
visible and their width increasingly narrows which is re- 
lated to the narrowing of the central peak as depicted 
in the inset. Obviously the central peak width is given 
by the low temperature fixed point value T c g which is 
already exponentially suppressed for X p h/F = 3. 

We have used a very large number of N z — 512 values 
in combination with a very small broadening b — 0.03 for 
the z-averaging of the spectral functions to ensure that 
the spectra depicted in Fig.[5]are really nearly broadening 
independent even at higher frequencies. 

The spectral peak at u> — remains pinned at its uni- 
versal value for T = in accordance with the Friedel sum- 
rule j^i^ 8 .!^ The spectral weight of this zero-frequency 
peak, however, is reduced to r o ff/To and redistributed 
to the higher energy phonon replica peaks. 

For X p h /To = 3, we checked in several very long NRG 
runs averaging up to N z — 1024 bath discretizations and 
using a very small broadening parameter of b = 0.01 
that the line width of the high energy phonon peaks are 
independent of the NRG broadening selected in Fig. [5] 

The increase of the spectral function and the maxima 
of the envelop functions at lo ±8Fo is related to the 
finite U e f{. Considering only the two local orbitals, d 
and c , of the IRLM J32]) which defines the initial NRG 
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FIG. 6. Feynman diagram of the electron-phonon self-energy 
E p (z) in Xp h , including the two attached incoming and out- 
going electron-propagator lines. The full line represents the 
local Green function Gd{z), the wiggled line the phonon- 
propagator. 




Q)/r n 



Hamiltonian Hq RG of the IRLM model, one can cal- 
culate the spectral functions exactly and finds single- 
particle excitations at wjv ~ ±U(N)/2. Since the value 
of U(N) can strongly depend on the iteration and ap- 
proaches its fixed point value U c s only for \uj\ <C r off , 
the peak positions correspond to the value of the renor- 
malized U(N)/2 where the NRG iteration N corresponds 
to the energy scale us ~ ujn oc A - ^^ 1 '/ 2 . This is in 
complete analogy to the single-impurity Anderson model 
where the bare high-energy value of U defines the loca- 
tion of the high-energy charge excitation peaks, while the 
low-energy fixed point value of C/ S °f m ~ Tr as shown by 
Hcwson et aiM- 

Although Fig. IDJa) indicates that for X p h/^o = 3, 
U e g is significantly lower than for the smaller values 
Ap/i/ro = 1-5,2, no additional charge peaks are found 
for the latter parameters. This is related to the fact that 
the self-energy contribution E ph (z) is proportional to X 2 h 
in leading order. If the magnitude of the self-energy is too 
small, the renormalized resonant level width T(N) dom- 
inates the envelop function. Once g > 1, Tq evidently 
becomes significantly suppressed, and spectral weight is 
increasingly shifted to higher frequencies. The broad 
phonon side peaks become observable once T e g <C To 
and the central resonance lost much of its weight. 



2. Weak coupling regime 

In the weak coupling regime, g 2 <C 1, the spectral func- 
tion is only slightly altered by the electron-phonon cou- 
pling. Starting from the free Green function of the RLM, 
the second-order contribution to the self-energy is given 
by the Feynman diagram shownin Fig. El Evaluating this 
leading order text-book diagram™ yields the analytical 
expression 

3mS ph (w - iS) = n\p h [g(u) ) + /(w + wo)]pd(w + cj ) 

+ 7rA p/j5( w o) + /(^o - - uj ) 

(47) 



FIG. 7. (Color online.) Imaginary part of self-energy contri- 
bution Tj ph (uj — i5) for Ed — 0, Wo = Fo and various values 
of X p h in the weak coupling regime \ p h/Fo < 1 and positive 
frequencies. The dashed line depicts the analytic self-energy 
given by Eq. Q47p an d evaluated for the same temperature 
T = 8.5 x l(T 3 ro as the NRG results. NRG parameter: 
iV s = 1500, A = 1.6, N b = 40, N z = 8, b = 0.1. 



where f(ui) (g(ui)) is the Fermi (Bose) function. The self- 
energy essentially vanished for \uj\ < u>q and T — ¥ and 
acquires a significant value only for > luq. The sharp- 
ness of the increase is governed by the Fermi function: an 
energy quantum of at least one bosonic excitation energy 
wo is needed to open up the additional decay channel. 

In a conserving approximation, one would replace the 
phonon propagator by its fully dressed version which in- 
cludes the shift of the phonon frequencies as well as a 
finite life-time broadening due to the local particle-hole 
excitations. This effect of the charge susceptibility has 
been neglected in Eq. (|4T|) . If we include this effect we 
would need to convolute Qm'E ph (uj — i5) with the true 
phonon spectral function which has been replaced by 5- 
functions in the derivation of Eq. (|4T|) . Therefore, we 
expect (i) a broadening of the steep increase above uiq 
with increasing X p h as well as additional features at mul- 
tiples of the phonon frequency, since S ph (w — iS) will also 
enter self-consistently the charge susceptibility. 

In Fig. the self-energy contribution 3ml] ph (u; — 
iS) close to the real axis is depicted for few values 
(Aph/To) 2 1. Dividing out the leading order pref- 
actor, A 2 ^, the magnitude of self-energy becomes nearly 
independent of X p h ■ The analytical expression of Eq. (|4"T)) 
added as dashed line shows a remarkably good agreement 
with the NRG self-energy obtained via Eq. (|45|) . Due to 
the broadening, the NRG self-energy cannot generate the 
same steep increase of the self-energy at u> « ojq. Note, 
that a shoulder is slowly developing at lo ~ 2ujq with in- 
creasing Xph, a precursor of additional phonon-peaks in 
the spectra. This is due to higher order processes not in- 
cluded in the second-order perturbation expansion, but 
present in a proper conserving approximation for S ph (z). 
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FIG. 8. (Color online.) Equilibrium spectral function for 
the molecular impurity Hamiltonian H m , Eq. (0 from weak 
to strong coupling in the anti-adiabatic regime too = 2Tq. The 
solid lines in (b) show spectra at T = 0.2Fo, the dashed lines 
in (a) the spectra for T = 1.2Fo. With e = and no = 0, E d 
is given by Ed = —g 2 Wo in Hi mp . NRG parameter: as in Fig. 

m 



FIG. 9. (Color online.) Comparison between the equilibrium 
and non-equilibrium spectral function for a single lead and 
the phonon energy ujq = 2To. The dotted lines show the 
equilibrium spectral function, the solid lines with the same 
color the spectra obtained for times t — > oo after a quench 
by switching on the electron-phonon coupling at time t = 0. 
NRG parameter: iV s = 1500, A = 1.6, N b = 40, N z = 8, 
6 = 0.1, T — 8.5 x 10~ 3 r . 



C. Particle- hole asymmetry 

Recently, Hiitzen et al.— investigated the I-V charac- 
teristics of the spin-less Anderson-Holstein model using 
an iterative real-time path integral approach,- 4 Within 
this approach the current is directly calculated using a 
generating functional, and no spectral functions are re- 
quired. Signatures of a Franck-Condon blockade have 
been reported using the local molecular Hamiltonian H m 

H m = oj^b + \ ph (tf + b)d)d + eSd 

defined in Eq. (JT]) for the parameters e = and uio — 
2r : the suppression of the current 1 ^ has been found 
when increasing the electron-phonon coupling X p h . Since 
Himp is identical to H m by setting uq = and Ed — 
e — E p , the suppression of the current at zero bias can 
be partially understood in terms of the polaronic shift of 
the single-particle energy Ed = e — g 2 uj - Increasing \ p h 
decreases Ed for constant e and, therefore, increases the 
particle-hole asymmetry. As a consequence, the stronger 
the electron-phonon coupling, the more energy is required 
to depopulate the level. 

The spectral functions are shown for two different tem- 
peratures and e = and ujq = 2Tq in Fig. [51 With 
increasing X p h, the main resonance peak is shifted to 
lower energies by E p . Additional phonon-side peaks oc- 
cur asymmetrically around the main resonance at Ed- 
The decrease of the spectra at uj = is clearly visible 
translating immediately to the reported^ suppression of 
the zero bias conductance. By comparing Fig. [HJa) and 
(b), we also note an increasing shift of the spectra to- 



ward the chemical potential and a visible narrowing of 
the main resonance when lowering the temperature for 
fixed X p h- 

D. Benchmark for non-equilibrium spectral 
functions 

The scattering states NRG approach to quantum 
transporter— relies on the analytically known density 
operator of a non-interacting model at finite bias4£ It 
is evolved to the density-matrix of the fully interacting 
problem using the TD-NRG i 40 i 41 This non-equilibrium 
density operator for the limit t — > oo is used to calcu- 
late the finite-bias retarded spectral function entering the 
equation (|24[) for the current-voltage characteristic of the 
junction^ 

The restriction to the single-lead version of the Hamil- 
tonian defined in Eq. |T]) in the absence of a bias can 
be used to benchmark the quality of the non-equilibrium 
spectra function algorithm and its limits. Starting from 
a decoupled electron-phonon system, \ p h — 0, and let- 
ting the system evolve using the full Anderson-Holstein 
model H = Hq + H p h, we must recover the equilibrium 
dynamics of H for t — > oo. Therefore, we compare the re- 
sults for the spectral function obtained from the standard 
equilibrium NRG algorithm 54 with the spectra obtained 
from a extension of the TD-NRG to spectral functions. 43 

The results for wo = 2Fo , Ed = and various electron- 
phonon coupling constants X p h are depicted in Fig. O 
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The dotted line shows the equilibrium spectral function, 
the solid line the spectra obtained from a quench by 
switching on the electron-phonon coupling at t = 0. The 
agreement remains very good even at high energies up to 
moderate values of X p h . 

The larger the phonon frequency, the better the agree- 
ment between the direct calculation of equilibrium spec- 
tra and the spectra obtained by the non-equilibrium ap- 
proach. For small ujo <Tq, however, the non-equilibrium 
approach has its limitations. The number of bosonic Fock 
states which are required even in equilibrium grows sig- 
nificantly. For our equilibrium calculation, we needed 
already Nb = 400 states, in a recent study— the rather 
very generous choice of NRG parameter, iVj, = 3000 and 
N a = 8000, have been used in a single lead model but no 
z-averaging was required for the thermodynamical prop- 
erties investigated in that paper. For the required z- 
averaging over 64 discretizations in an TD-NRG calcula- 
tion, this would be computationally much to expensive 
and is, therefore, beyond the reach of our approach. 

V. QUANTUM TRANSPORT: RESULTS 
A. Scattering-states NRG 

This section is devoted to the results of the non-linear 
I-V characteristics calculated using the SNRG. The dif- 
ferential conductance G(V) — dl/dV has been obtained 
numerically from the I(V) curves. 

Section IV A 11 is devoted to the particle- hole symmet- 
ric regime of the model, i. e. Ed = 0. We investigate the 
evolution of the I-V characteristics from a junction with 
symmetric couplings (R = 1) to the tunneling regime 
(R — s- oo.) In Sec. IV A"2l we focus on the PH-asymmetric 
model and investigate the Franck-Condon blockade for 
e = Ed — E p = and a finite electron-phonon coupling. 

1, Particle-hole symmetric model 

The electrical current calculated via Eq. ([24]) is shown 
as a function of the applied bias in Fig. ITUT a) for three 
different ratios of R — Tl/Yr, a phonon frequency uiq = 
2Tq and a fixed electron-phonon coupling g = 1. For 
better comparison, the prefactor Go defined in Eq. (I25[) . 
has been divided out. Go only contains the trivial reduc- 
tion of the current upon increasing of R starting from 
its maximum of e 2 /h for R = 1. For a PH-symmetric 
Hamiltonian, the current remains PH-symmetric even for 
an asymmetric junction as can be seen analytically from 
Eq. (21. 

The differential conductance G(V) = dl/dV depicted 
in Fig. [TUTb) has been obtained by numerically differenti- 
ating the data of Fig. HOT a) . For a better comparison, we 
added the corresponding equilibrium transmission func- 
tion obtained from the spectral function shown in the 

Fig.m 
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FIG. 10. (Color online.) (a) I(V) for the symmetric model, i. 
e. Ed = 0, and three different ratios R = Tl/Tr — 1, 10, 100. 
The vibrational parameters are luq = \ p h = 2ro. (b) dl/dV 
in units of Go as function of eV obtained by numerically dif- 
ferentiating the curves in (a). Additionally, the equilibrium 
transmission function T(uj) = nTopd{w) obtained from the 
spectral function pd(oj) of Fig. [9] has been added for compar- 
ison. NRG parameter: N s = 2000, A = 2, N b = 35, N z = 64, 
6 = 0.15, T = 6.5 x 10" 3 r . 



Particle-hole symmetry and a Fermi-liquid equilibrium 
ground state yields a pinned peak at zero bias which ap- 
proaches Go for T — > which can be understand quite 
easily by applying the Friedel sum rule i 38 ' 39 The narrow- 
ing of the central resonance due to the renormalization 
of To — > r o ff, exemplified in Fig. [S] prevails also in the 
differential conductance. 

Qualitatively G(V) traces the transmission function 
T(uj) = ttTqp(uj) but quantitatively significant differ- 
ences are observed. The additional transmission max- 
ima are clearly visible at finite voltage which are related 
to phonon-assisted tunneling increasing significantly the 
current above a threshold. This threshold position de- 
pend on the asymmetry R of the junction. 

We expect that G^dl /dV approaches the equilibrium 
transmission function in the tunneling regime R — > oo 
and T — > 0. This is clearly the case for R = 100 and 
low voltages. We note, however, that the peak position 
of first maxima at eV ~ ±wo is well reproduced but the 
peak height is smaller than expected from the equilibrium 
transmission function T{uS) = 7iTo/9,i(u;). This might be 
caused by remaining non-equilibrium effects still relevant 
for R = 100. However, we believe that this is cause by 
the limitations of the numerical accuracy of the SNRG at 
large bias caused by the discretisation and the broadening 
of the spectral function in combination with the TD-NRG 
time evolution. 

By decreasing the asymmetry from R = 100 to R = 10 
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FIG. 11. (Color online.) (a) I(V) for the particle-hole 
asymmetric model, Ed /To = — 2(e = 0) and an asymmet- 
ric junction 7? = Tl/Tr = 10. The vibrational parameter are 
ujo = 2ro, and \ v h = 2Fo. The solid line shows the SNRG cur- 
rent, The black dotted line indicates the \ p h = £ = current, 
the blue dashed line current for a shifted level E 3 /Tq — —2.5 
in an non-interacting junction (RLM)at T — 0. (b) dl/dV in 
units of Go as function of eV obtained by numerically differ- 
entiating the SNRG curve in (a). NRG parameter: as in Fig. 



small differences are visible which can be traced back 
mainly to the changes chemical potentials [i L and 
for the same voltage drop. For a symmetric junction 
R = 1, the strongest non-equilibrium effects are ob- 
served. G(V) does not follow the expectation G(V) oc 
[T(eV/2) + T(— eV/2)] where the equilibrium transmis- 
sion function T(uj) naively has replaced full bias depen- 
dent spectra pd{u, V) in Eq. (j24l) . The width of the cen- 
tral peak is smaller that expected for G(V) oc [T(eV/2)] 
and the peaks from the phonon assisted tunneling occur 
closer to ±2wo as suggested by the equilibrium trans- 
mission function. This is due to a significant change 
of the non-equilibrium spectral function with increasing 
bias voltage. Spectral weight is redistributed to higher 
frequency due to backscattering processes for which ad- 
ditional phase space becomes available at finite voltage. 

The full width at half maximum (FWHM) of the zero 
bias dl/dV peak is not simply given by 2r e ff as sug- 
gested by a simple generalization of equilibrium spectral 
function. Non-equilibrium effects cause a bias depen- 
dency of the spectral function: the width of dl/dV is 
significantly smaller than predicted from such a naive 
[T(eV/2) + T(-eV/2)] fit to G(V) for R = l. 



2. Particle-hole asymmetric model 

For a particle-hole asymmetric regime of the model 
and asymmetric lead couplings, the current as function 



FIG. 12. (Color online.) (a) I(V) for E d /T = -0.6(e/T = 
1.65) and R = Tl/Tr — 10. The vibrational parameter are 
cjo = 4Fo, and X p h = 3Tq. The solid line shows the SNRG 
current, The black dotted line indicates the X p h = £ = 
current (RLM) at T = 0. (b) dl /dV in units of Go as function 
of eV obtained by numerically differentiating the SNRG curve 
in (a). NRG parameter: as in Fig. 1101 



of the applied bias voltage and corresponding the dif- 
ferential conductance G(V) is depicted in Fig. QT] We 
have set e = in H m in resonance with both chem- 
ical potentials at zero bias, the phonon frequency to 
wo/ro = 2 and selected a moderate electron-phonon cou- 
pling g = Aph/uiQ = 1. As discussed before, the bare 
single-particle level e is shifted by the polaron energy 
to = —E p = — 2T . The corresponding equilibrium 
spectral function has been plotted in Fig. |S]which clearly 
documents the significant asymmetry. 

The strong suppression of the current with increas- 
ing Xph and fixed e, as seen by comparing the X p h = 
and Xph /To = 2 curves has been interpreted as Franck- 
Condon blockade physics. 19 Clearly, the leading order ef- 
fect stems from a polaronic shift E p of the single-particle 
level. The SNRG curve tracks the current through a 
shifted level at Ed/To = —2.5 and X p h — rather well 
for positive voltages as indicated by the read dashed line. 

We have already seen in the equilibrium spectral func- 
tion - green curve in Fig. [5] - that the electron-phonon 
interaction causes a redistribution of spectral weight be- 
low the chemical potential. Each phonon peak contains 
only of a fraction of the spectral weight, and spectral 
weight has been shifted even below — 10Fo which con- 
tribute only at much larger negative voltages. In the 
asymmetric junction (R — 10), the current \I\ increases 
significantly for eV < —E p , However, the magnitude of 
current clearly remains lower than the reference current 
(dashed line) for a non-interacting shifted level. This is 
consistent with the the qualitativ features of the equi- 
librium spectra and leads to the Franck-Condon current 
suppression of the current at small and intermediate bias 
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voltage. A second phonon induced maximum in G{V) is 
found below V < —E p —ujq and corresponds to the second 
peak in the spectral function. A very shallow maximum 
at eV/To » — 7 could be related to the third phonon peak 
in the spectral function. We believe, that the small oscil- 
lations in die dl/dV curve for positive bias voltage has 
no physical significance and is caused by the numerical 
differentiation of the numerical I-V data. 

Maintaining the same coupling asymmetry R = 10 but 
increasing the phonon energy to loq/Tq = 4, the current 
vs voltage is depicted in Fig. IT27 a) and the corresponding 
dl/dV curve in Fig. [T27 b) for a moderate electron phonon 
coupling g — 0.75. We also added the equilibrium spec- 
tral function to the (b) as a comparison. Since R = 10, 
dp(uj,V) /dV remains significant, and dl/dV does not 
trace the equilibrium spectra p(uj, V = 0). There remains 
a significant renormalization of the spectral function for 
a moderate asymmetry which only must vanish in the 
limit R — y oo or R — y 0. 



B. Tunneling limit for the crossover regime 

In the adiabatic regime the physics is dominated by the 
charge-fluctuation scale Tq being the largest local energy 
scale in the problem. The weak electron-phonon coupling 
causes only small renormalization of the electronic and 
phononic degrees of freedom. In this regime, Keldysh ap- 
proaches have been successfully applied to the problem: 
the phonon frequency is reduced by particle-hole exci- 
tations in leading order which is well captured already 
in conserving one-loop diagrams^^i for the electron and 
the phonon propagator. 

Recently, it was shown^ 2 - that the anti-adiabatic regime 
extends from luq > Tq to luq w 0(T), as long as 
the polaronic shift E p exceeds the charge-fluctuation 
scale. The crossover regime from this extended anti- 
adiabatic regime to the adiabatic regime where the po- 
laron contains less and less phononic excitations has been 
defined 2 ^ by the parameter hierarchy E p > Tq > luq and 

r c ff ~ w . 

Although this regime is accessible to the equilibrium 
NRG, a very large number of phononic states are needed 
for accurately tracking the equilibrium flow/22, The SNRG 
relies on the switching on the electron-phonon coupling 
and let the system evolve from a X p h = to a finite \ p h- 
Apparently, the TD-NRG is not able to reproduce reli- 
ably the equilibrium Green functions for larger electron- 
phonon couplings in the crossover regime. This might 
also be related to the logarithmic discretisation of the 
scattering states in combination with the change of the 
ground state s 42 ' 56 

Therefore, we have focused on the tunneling regime 
(R — > oo) where the retarded Green function becomes 
bias independent and V) can be replaced in Eq. (|24|) 
by pd(uj,V — 0). We have included Nf, — 300 phonon 
Fock states and kept N s = 2500 NRG after each iteration 
to calculate the equilibrium spectra using only a single- 




FIG. 13. (a) I(V) for the particle-hole asymmetric model, 
Ed/To = 1.2 in the crossover regime. The phonon energy 
loo = 0.26r , and A ph = 0.6IV (b) Ad 2 //dV 2 in units of Go 
as a function of eV obtained by taking analytically the second 
derivative of Eq. (|24[) in the limit R — ¥ oo. NRG parameter: 
N s = 2500, A = 1.8, A^ b = 300, N z = 8,6 = 0.1, T = 8 x 

io- 4 r„. 



lead. 

Since the charge-fluctuation scale r is the dominat- 
ing energy scale, G(V) is unsuitable to reveal subtle 
inelastic processes induced by the electron-phonon cou- 
pling in the cross-over regime. The information about 
the inelastic processes, however, is encoded in the second 
derivatively of the I(V) curve, d 2 I/dV 2 . 

For the tunneling I(V) is depicted in Fig. [ToT a). Al- 
though, luq < To, the polaron shift is given by E p /Tq = 
1.38 and exceeds Tq, while the renormalized r o ff remains 
above loq. Therefore, the choice of parameters lies in 
the crossover regime between the extended anti-adiabatic 
and the adiabatic regime. 

Since d 2 //dV 2 is up to some smaller voltage depen- 
dence proportional to the derivative of the spectral func- 
tion, it is sensitive to the derivative to the phonon self- 
energy. Although, this absolute magnitude is small, its 
derivative can be huge if it contains a threshold set 
by a Fermi-function shifted by one or multiples of the 
phonon frequencies. In second order Keldysh approxima- 
tion, these steep features dominate d 2 I /dV 2 on voltages 
eV/ujQ « 0(1) as seen in Fig. 11 of Ref. 10. The sign 
change of the d 2 //dy 2 for small voltages can only occur 
for very pronounce changes in the self- energies. 

In the NRG, the features are much less dominant: the 
method has its limitation due to the discretization and 
the artificial broadening which puts limit to the gradient 
of any self-energy change as already discussed above. In 
order to extract some information about inelastic scat- 
tering processes, we have calculated d 2 /o/dV r2 by using 
a Lorentz fit to the spectral function centered around 
Ed which is the mean-field (MF) reference value. Now 
we calculate the full d 2 I/dV 2 and define Ad 2 I/dV 2 = 
d 2 I/dV 2 — d 2 /o/dT^ 2 as difference between the full sec- 
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ond derivative and the MF result. This is plotted in 
Fig. nUT b). Now we notice the sharp features around 
eV = ±wo and at multiples of the phonon frequencies. 
This is clearly visible in the inset of Fig. [TSTb) where 
the bias is measured in units of ujq- While for X p h — > 0, 
the self-energy is dominated by threshold of the inelas- 
tic processes at |w| = ujq, for increasing X p h multiple 
phonon absorption and emission processes start to con- 
tribute to the self-energy. Similar to Fig. [7] where two 
phonon processes are clearly visible, this multi-phonon 
processes at moderate coupling yield the additional fea- 
tures in Ad 2 I/dV 2 . The one and two-phonon processes 
have been also seen in a much more pronounce manner 
in the literature^ using Keldysh perturbation theory. 

VI. DISCUSSION AND OUTLOOK 

We have extended the scattering states NRG to the 
charge-transport through a molecular junction. To set 
the stage for the non-linear transport results, we have 
provided a detailed analysis of the low-temperature equi- 
librium physics of the spin-less Anderson-Holstein model 
from an NRG perspective. We have shown that the 
model can been mapped onto an interacting resonant 
level model in the extended anti-adiabatic regime and 
extracted the renormalized charge-transfer scale T e g and 
the effective Coulomb interaction U e ff in a regime com- 
plementary to the recent study of Eidelstein et al.. 22 

We have calculated the equilibrium spectral func- 
tions and used those to benchmark our non-equilibrium 
algorithm^ We have demonstrated that for the weak 
coupling limit, (X p h/ujo) 2 *C 1, the NRG tracks perfectly 
the self-energy obtained from the second-order Feynman 



diagram. For large electron-phonon couplings the typi- 
cal phonon replica peaks of the exact atomic solutio n 15 i 28 
are found. Using extensive z-averging the spectra become 
independent of the NRG broadening even at higher fre- 
quencies. 

Exemplified in Fig. [TU1 the reduction of the charge- 
transfer scale r o ff due to the electron-phonon coupling 
in the equilibrium properties conveys into a narrowing 
of the differential conductance. In contrary to a recent 
QMC studyi^ for the spin-full model with Coulomb inter- 
action, our results clearly show the phonon side peaks ex- 
pected from the equilibrium spectra. The location, how- 
ever, depends on the two chemical potentials as well as 
the coupling asymmetry of the junction. 

Gating the junction away from the particle-hole sym- 
metric point reveals the Franck-Condon blockade physics 
with increasing electron-phonon coupling. The leading 
order effect is understood in terms of a polaronic level 
shift, and the I(V) curve tracks a shifted resonant level 
model for positive bias voltages. For negative bias, how- 
ever, the current is suppressed due to a redistribution of 
spectral weight to higher frequencies. 
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